Towards an understanding of the stability properties of 
the 3+1 evolution equations in general relativity 



O 

o 
o 

(N 
-i— > 
O 

O 

(N 
> 
On 

o 

00 

O 

On 
ON 

"o 
cr : 

u '. 



Miguel Alcubierre^, Gabricllc Allen^, Bernd BriigmannW, Edward Scidcl^ 1 ' 2 ), and Wai-Mo Suen^ 3 ' 4 ) 

Max-Planck-Institut fur Gravitationsphysik, Albert- Einstein- Institut, Am Muhlenberg 1, 14476 Golm, Germany 
' 2 ' National Center for Supercomputing Applications, Beckman Institute, 405 N. Mathews Ave., Urbana, IL 61801 
^ Department of Physics, Washington University, St. Louis, MO 63130 
Physics Department, Chinese University of Hong Kong, Hong Kong 
(February 4, 2008; AEI-1999-19) 

We study the stability properties of the standard ADM formulation of the 3+1 evolution equations 
of general relativity through linear perturbations of flat spacetime. We focus attention on modes 
with zero speed of propagation and conjecture that they are responsible for instabilities encountered 
in numerical evolutions of the ADM formulation. These zero speed modes are of two kinds: pure 
gauge modes and constraint violating modes. We show how the decoupling of the gauge by a 
conformal rescaling can eliminate the problem with the gauge modes. The zero speed constraint 
violating modes can be dealt with by using the momentum constraints to give them a finite speed 
of propagation. This analysis sheds some light on the question of why some recent reformulations of 
the 3+1 evolution equations have better stability properties than the standard ADM formulation. 
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I. INTRODUCTION 

There has been an intense effort in trying to de- 
velop Numerical Relativity for the study of astrophysi- 
cal phenomena involving black holes and neutron stars. 
Most investigations in numerical relativity for the last 
30 years have been based on the Arnowitt-Dcscr-Misncr 
(ADM) yj] system of evolution equations and many im- 
portant results have been obtained in spherical symme- 
try and axisymmetry. However, in the general three- 
dimensional (3D) case which is needed for the simula- 
tion of realistic astrophysical systems, it has not been 
possible to obtain long term stable and accurate evolu- 
tions (although some good progress has been made, see, 
e.g., ||-|^1). One might argue that present day compu- 
tational resources are still insufficient to carry out high 
enough resolution 3D simulations. However, the diffi- 
culty is likely to be more fundamental than that. There is 
no theorem guaranteeing the well-posedness of the initial- 
boundary value problem for the full ADM system. In par- 
ticular, one must consider the possibility that free evo- 
lutions using the ADM system might be unstable, e.g., 
against constraint violations in 3D. There are also well- 
known complications due to the gauge (coordinate) de- 
grees of freedom in the theory. This is one of the major 
open problems in numerical relativity. 

Against this background of need and failure to obtain 
long term stable, accurate 3D simulations in numerical 
relativity, in the last decade there has been a lot of effort 
looking for alternative formulations of the 3+1 equations, 
which can be roughly separated in two directions. (I) 
In the mathematical direction, several first order hyper- 
bolic formulations have been proposed, and conditions 



on well-posedness of the initial-boundary value problem 
have been 

studied [|-23|. Unfortunately there is as yet no evi- 
dence that the hyperbolic re-formulations lead to signifi- 
cant improvements in general 3D numerical calculations 
(despite encouraging results in the spherical symmetric 
case 

@)- (II) In the more "empirical" direction there 
have also been various attempts to improve stability and 
accuracy by modifying the ADM system. To avoid in- 
stabilities due to constraint violation, fully or partially 
constrained evolutions have been tried and the addition 
of "constraint enforcing terms" into the ADM evolu- 
tion equations has been proposed and attempted |pi 25 
(cmp. j26|). Methods to better enforce gauge conditions 
have also been suggested J2?| . Most significant and rele- 
vant for our present paper is an approach based on sepa- 
rating out the conformal and traceless part of the ADM 
system, first developed by Shibata and Nakamura [EsJ . 
Unfortunately, the strength of the Shibata-Nakamura ap- 
proach was not widely appreciated, until Baumgarte and 
Shapiro [^9| compared the standard ADM formulation 
with a modified version of the Shibata-Nakamura for- 
mulation on a series of test cases, showing the remark- 
able stability properties of the conformal-traceless (CT) 
system. This has triggered much recent research in the 
community, including what we are reporting here and 
in a companion paper. There also have been interesting 
results connecting the conformal approach to the hyper- 
bolic approach @|o|j3lJ. 

In this and a companion paper we compare the stan- 
dard ADM equations to the CT equations of Shibata- 
Nakamura and Baumgarte-Shapiro in different imple- 
mentations. In the companion paper p2j, we show em- 
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pirically the strength of this system over the standard 
ADM equations in numerical evolutions, at least in some 
of the implementations of the former set of equations. 
We study in particular the CT formulations in numerical 



evolutions of strongly gravitating systems (see also |33|) 
and when coupled to hydrodynamic evolution equations, 
extending previous studies of weak fields |^9| and of pre- 
determined hydrodynamic sources |34j . The main conclu- 
sion is that the CT formulation is more stable than the 
standard ADM formulation in all cases, while it needs 
more resolution for a given accuracy than ADM in some 
cases. 

In this paper, we aim at developing a mathematical un- 
derstanding of the stability properties of the ADM and 
the CT equations. Ideally one would like to know if the 
different systems are well-posed. However, the systems 
of equations as they stand are mixed first-second order 
systems and as such are not hyperbolic in any immediate 
sense. This makes a study of their well-posedness par- 
ticularly difficult. Because of this fact, we have chosen 
instead to study linear perturbations of a flat background 
and do a Fourier analysis. We believe that this analysis, 
though only valid in the linear regime, reveals important 
information about the stability properties of the different 
formulations. 

We study in particular two types of zero speed modes 
that appear in the standard ADM formulation, the 
"gauge modes" and the "constraint violating modes" , 
and what they turn into in different implementations of 
the CT system. The main result of this paper is a con- 
jecture that the zero speed modes are responsible for the 
instabilities seen in the integration of the ADM system, 
and a suggestion of how they could be handled to ob- 
tain stable evolutions. We stress the point that we do 
not believe that these instabilities are of numerical ori- 
gin. Instead, we believe that they correspond to genuine 
solutions of the exact system of differential equations. A 
related analysis to the one we present here, but along 
different lines, was recently carried out by Frittelli and 
Reula ff|. 

In sectio n O , we study the linearized ADM equations. 
In section |l| we introduce a model problem to help us 
understand the effect of the zero speed modes. In sec- 
tion |^ we discuss the gauge modes, and in section |v| 



numerical 

We conclude with section fvTJ 



the constraint violating modes. In section VI 
examples are considered 
Comments on finite difference approximations to the lin- 
earized ADM equations can be found in the appendix. 

A final comment about the language used to describe 
the solutions to the different systems of equations. We 
have chosen to refer to all solutions that satisfy the con- 
straints as physical solutions, and those that do not as 
unphysical. According to this criterion we will consider 
pure gauge solutions as physical solutions, even if they 
contain no real physical information. 



II. THE LINEARIZED ADM EQUATIONS 

Let us consider first the standard ADM evolution equa- 
tions for the spatial metric and extrinsic curvature K%j 
which in vacuum take the form: 

(d t - Cp) Qij = -2aKij, (I) 
(dt - Cp) Kij = —DiDja + a (Rij + KK l3 (2) 



with Cp the Lie derivative with respect to the shift vector 
[3 l , a the lapse function, Di the covariant derivative with 
respect to the spatial metric, Rij the Ricci tensor of the 
3-geometry, and K — g lJ K^. 

Together with the evolution equations, one must also 
consider the Hamiltonian constraint 



R + K 2 - K i0 K %i =0, 
and the momentum constraints 

Dj (K ij — g ij K) = 0. 



(3) 



(4) 



Let us now take geodesic slicing a — 1 and zero shift 
[3 l = 0, and consider as well a linear perturbation of flat 
space 



6>ij -\- h'ij , 



(5) 



with hij « 1. The evolution equations now reduce to 



d t hi 



-2 K, 



d t K ij =R%\ 



(6) 
(7) 



where the linearized Ricci tensor is given by 

R$ = -1/2 (vLfctf - diTj - djTi) . (8) 
and where we have defined (h = tr/iy) 

Ti :=J2 d khik- 1/2 dih, (9) 

k 

Notice that I\ is just the linearized version of g mri T l mn . 

In the same way, we find that the linearized approxi- 
mation to the constraints is 



^""^ dkfk — 0, (hamiltonian) 



where now 



dtfi — 0. (momentum) 



/, := ^dkhik - dih. 



(10) 
(11) 

(12) 



The structure of the constraints is quite interesting. 
They just state that the vector / should be both diver- 
genceless, and time independent. Notice that both these 
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conditions would be trivially satisfied if we were to choose 
/=0, which somewhat counter-intuitively amounts to 
three conditions instead of four. Notice also that asking 
for a transverse (dkhik=0) and traceless (h=0) solution 
means that f—0, so the constraints are satisfied automat- 
ically. This is precisely what is done when one chooses 
the standard transverse-traceless (TT) gauge. 

Having found the linearized evolution equations, we 
now proceed to do a Fourier analysis. Without loss of 
generality (but see Appendix), we can take the plane 
waves to be moving in the x direction. The result for 
any other direction can be recovered by a simple tensor 
rotation later. We then assume that we have a solution 
of the form 



a i{uit—kx) 



Equation (^) implies 

= - (iuj/2) hij. 
Substituting this into Eq. (Q) we find 
uj 2 h. = k 2 Mh, 
where we have defined the six-dimensional vector 



h := (h xx , 
and the matrix 



foyy-! hzzj h X y 7 h X z? ^lyz 



(13) 
(14) 

(15) 

(16) 

(17) 



M 



/ 1 1 \ 

1 

1 





\ 1 ) 

One can also find that the constraints ( |To| ) and (|l 
reduce to the three (not four!) conditions 



(18) 



b yy 



h zz =Q, 



h xz = , 



(19) 
(20) 
(21) 



where the first one of these equations results from both 
the hamiltonian constraint and the x component of the 
momentum constraint, and the last two result from the 
y and z components of the momentum constraint respec- 
tively. 

It is straightforward to calculate the eigenvalues A and 
eigenvectors of the matrix M. They turn out to be 



• A = 0, with corresponding eigenvectors 
V! = (1,0,0,0,0,0), 



V2 

v 3 



(0,0,0,1,0,0), 
(0,0,0,0,1,0). 



(22) 
(23) 
(24) 



• A = 1, with corresponding eigenvectors 

v 4 = (2, 1,1, 0,0,0), (25) 

w B = (0,1,-1,0,0,0), (26) 

v 6 = (0,0, 0,0, 0,1). (27) 



What sort of solutions do the different eigenvectors 
represent? There are four different types of solutions: 

1. Two modes that satisfy all the constraints that 
travel with the speed of light (A = 1) represented 
by the transverse-traceless vectors v§ and Vq. 

2. One mode that violates both the hamiltonian and 
the x component of the momentum constraints 
(compare with Eq. (19)), that also travels with the 
speed of light (A = 1) represented by the vector V4. 

3. Two transverse modes that violate only the mo- 
mentum constraints (compare with Eqs. j20| ) 
and j2l|)), and "travel" with speed zero (A = 0) 
represented by the vectors V2 and V3. 

4. One mode that satisfies all the constraints that also 
has speed zero (A — 0) represented by the vector 

Vi. 

The three constraint satisfying modes are clearly phys- 
ical solutions. Of these, the two transverse-traceless trav- 
eling modes (us and vq) correspond to the standard gravi- 
tational waves. What is the remaining physical mode v±? 
The only possibility is for it to be a pure gauge mode. To 
see that this is indeed true all we need to check is that it 
corresponds to a solution for which the 4-curvature Rie- 
mann tensor vanishes. For this we start from the Gauss- 
Codazzi relations, which to first order are 



(4) pm _ (3) r>m 

Hjk ijk ' 

^R% k = d k K t3 - djK lk 

(4)d0 —_S).V 
n i0k — U * /V ' 



(28) 
(29) 
(30) 



Now, the fact that v\ has dependence only on x (by 
construction), and has a component corresponding only 
to h xx implies that the r.h.s. of ( p9| ) vanishes and hence 
^R^ k —0. Also, since this mode has zero speed, it corre- 
sponds to a mode for which dtKi k —0 which in turn means 
that ^R^ ok =0. Finally, it is not difficult to see that for a 
solution that depends only on x and for which only h xx is 
non-zero, the 3-curvature vanishes as well, which tells us 
that ^R"j k =0- The 4-Riemann is then identically zero, 
so the mode v\ is a pure gauge mode. 

The presence of the zero speed modes (v±, V2 and V3) is 
troublesome: They do not represent non-evolving modes 
as one might think at first sight, but rather they rep- 
resent modes that annihilate the Ricci tensor. As such, 



3 



they correspond to solutions for which the extrinsic cur- 
vature remains constant in time, and the metric func- 
tions grow linearly (the linearly growing gauge modes 
have been studied before in |35|,|36|] ) . With the full non- 
linear ADM equations, this linear growth is likely to lead 
to an instability. 

In the next section we use a simple model problem to 
show how zero speed modes can indeed become unstable 
in the presence of non-linear terms. 



III. ZERO SPEED MODES: A MODEL PROBLEM 

To understand the effects of zero speed modes on sta- 
bility, we study the simple case of the one-dimensional 
wave equation with a non-linear source term F: 



Evolution of 



d t 2 <t>-ed x 2 4> = 5F{<t>,d t ci )l d t( t>). 



(31) 



We investigate the stability of the system for different 
values of e and <5. We will call the system unstable if the 
magnitude of <j) grows faster than exponential in time at 
any fixed value of x, and stable otherwise. 

For 6 = 0, but e not equal to 0, we have the usual wave 
equation. A Fourier decomposition of the form used in 
the last section reveals two eigen-modes with propagation 
speeds A = ±y/e. The system is stable for all values of e 
including zero, if there is no source term ((5 = 0). With 
a source term, it will still be stable for non-zero e, but 
not so if e becomes zero. For zero e, the two propagation 
speeds degenerate to zero, and the system is unstable for 
a general source term. This can be shown analytically by 
writing (|3l|) in first order form: 




where D 



0, 











-& 


(32) 








d t cj>. 


For e not 


equal 



and 7r 

to 0, the characteristic matrix (the matrix multiplying 
the d x term above) has three independent eigenvectors: 
(1, 0, 0), (0, 1, y/e), and (0,1,—^). The eigenvector ma- 
trix and its inverse have bounded norms. The system is 
therefore str ong ly hyperbolic, which in turn guarantees 



its stability |37[ . 

When e = 0, two of the eigenvectors become degen- 
erated and the system becomes weakly hyperbolic [ p7[ . 
For 6 = 0, the system is still stable, with at most linear 
growth in <p. But for 6 non-zero and with a general source 



term, the system is unstable |37||3b 

A2 



As an example, we take F = <\> in (|E 



d t 2 ci>-ed x 2 4> = 6 4> 2 - 



(33) 



When e is non-zero, there are no zero speed modes and 
the evolution is stable. In Fig. |l| we show the evolution 
of 4> at various times (from t = to t = 30 in equal time 




FIG. 1. Evolution of <j> described by Eq. (|3J), with e = 1 
and 8 = —0.01 at various times (from t = to t = 30 in equal 
time intervals). 



intervals) for the case of e = 1 and 6 = —0.01 (the initial 
data is a Gaussian wave packet). This evolution is very 
similar to that of a non-linear gravitational plane wave 
(see |§). 

Next we tune e down to zero in Eq. ([33]). The prop- 
agation speed of the eigenmodes becomes zero. The ini- 
tial Gaussian profile now does not propagate, instead it 
decreases in amplitude initially, becomes negative and 
eventually blows up. See Fig. || for the evolution up to 
t = 5, with the same initial data as before. In fact, this 
system is simple enough to be solved exactly. One can 
show that, at given value of x, the solution blows up as 
— \/{t — c{x)) 2 , where c is a constant depending on the 
initial value of </> at that point. From this it is clear that 
</> in fact becomes infinite after a finite time. 

We have studied examples with different source terms 
and have seen similar behavior, namely the systems be- 
come unstable when e goes to zero. To relate more di- 
rectly this scalar field instability to the ADM equations, 
we insert variable parameters (e's) into the linearized 
ADM system studied in the previous section. We exam- 
ine the case in which the matrix M (in Eq. (|l8|)) contains 
variable parameters t\, and €3: 

/eill 0\ 

1 

1 

e 2 

e 3 

\ 1/ 



M f 



(34) 



For non-zero (positive) e's, the corresponding set of 
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Evolution of 




FIG. 2. Evolution of 4> described by Eq. @, with e = 
and 8 = —0.01 at various times (from t — to t = 5 in equal 
time intervals). 



ferent from those in the zero-wave-speed case we focused 
on above. These "blowing up solutions" are blowing up in 
a global manner, and can come into existence in our nu- 
merical evolutions only if we choose boundary conditions 
that allow them. In numerical evolutions (at least those 
considered in this paper) we typically start the evolution 
at a certain initial time in a compact computational do- 
main with a certain chosen set of boundary conditions. 
The "blowing up solutions" , which are blowing up in a 
global manner can be excluded by an appropriate set of 
boundary conditions. On the other hand, in the case 
when e = and 6 > 0, the unstable solution involves an 
arbitrary function of x. One can see that any initial data 
with positive <p will cause a local blow up, independently 
of its initial profile. It cannot be excluded by choosing 
boundary conditions. The locality of the instability is 
the crux of the problem making it dangerous in numeri- 
cal evolutions. 

In the next sections, we focus on the zero speed modes 
in the case of the Einstein equations. We show how one 
can deal separately with the gauge mode and the con- 
straint violating modes. 



second order differential equations has no zero speed 
modes. To investigate its stability, it is straightforward 
to break this second order system into a first order system 
as in the scalar field study above. It is then easy to show 
that the resulting first order system is strongly hyperbolic 
and hence stable. As we turn the e's to zero (recovering 
the ADM system), zero speed modes appear in the sec- 
ond order system and the corresponding first order sys- 
tem becomes only weakly hyperbolic. This is precisely 
what happened in the scalar field example above. We 
hence conjecture that the existence of zero speed modes 
and the related weak hyperbolicity is at least one of the 
reasons why the ADM system becomes unstable in nu- 
merical evolutions when non-linear source terms cannot 
be neglected (i.e. unstable when gravity and/or gauge 
effects become strong). 

As an explicit example of such a transition between lin- 
ear growth and non-linear blow up in the ADM case, we 
note the well known case of focusing in geodesic slicing. 
It is precisely the zero speed gauge mode discussed in the 
previous section the one that represents the focusing of 
geodesic slicing. In the full non-linear case this focusing 
produces a coordinate singularity causing a blow up in a 
finite time. 

One last comment comparing different blowing up solu- 
tions is in order. We note that the non-linear wave equa- 
tion (|33j) described above has solutions that blow up in a 
finite time even in the case of a non-zero wave speed. For 



e — 1 and 8 — 1, two such solutions are: <fi = 
(with c a spatial constant) and 4> = —4/(t 2 — x' z ). How- 
ever, these "blowing up solutions" are fundamentally dif- 



-6/(t- 



IV. DEALING WITH THE GAUGE MODE: 
DECOUPLING K 

In trying to deal with the zero speed modes, we will 
first concentrate on the pure gauge mode: the mode V\, 
associated with h xx in the analysis of section |n]. Since 
this mode satisfies all the constraints, it represents a 
physical solution of the evolution equations (even if it 
only corresponds to a non-trivial evolution of the coor- 
dinate system), and hence cannot be eliminated. The 
most we can hope to achieve is to decouple it from the 
rest of the evolution equations, so that it will be immune 
to possible numerical errors, in particular those coming 
from the complicated Ricci tensor terms driving the evo- 
lution. 

Remarkably, such a decoupling is not difficult to 
achieve. Following [p9| , p8|j3C| l we first conformally rescale 
the metric in the following way 



9ij 



-4<f> 

e yij i 



(35) 



with (j) chosen in such a way that the rescaled metric c/ij 
has unit determinant, 



T2 l0g5 ' 



(36) 



We also define the conformally rescaled, trace- free part 



of the extrinsic curvature Ka as 



Aij — e 4 ^ ( Kij — —gijK 



(37) 



The ADM equations (|l|) and (|3|) can now be rewritten 
as the following system of 14 evolution equations 



5 



(d t - Ct 



6 



(d t - £p)gij = -2aAij, 
(d t -Cp)K= -D i D i a + a (R + K 2 ) 



(d t -C )Ai 



-4<p 



(—DiDja + aRi 



a ( KAa - 2A a A\ 



(38) 

(39) 
(40) 

(41) 

(42) 



subject to the extra constraints g = 1, trA = 0. 

The hamiltonian and momentum constraints now be- 



come as 



R - AijAv + 2K 2 /3 = 0, (43) 
Dj (Av - 2g^K/i\ + QA l ^dj(t> = 0. (44) 

Notice that now we have separated out the "gauge" 
variables {(f), K}, but we haven't yet decoupled the evo- 
lution equation for K from the Ricci tensor. This last 
step can be achieved by making use of the hamiltonian 
constraint above. Doing this, we can eliminate all refer- 
ence to the Ricci tensor from the evolution equation for 
K . One can also use the hamiltonian constraint to elimi- 
nate the Ricci scalar from the evolution equation for A^ . 
In fact, one can consider adding an arbitrary multiple of 
the hamiltonian constraint to this equation. We will then 
consider the evolution equations 



(d t -£?)<!> = -~K, 
(d t - Cpjgij = -2aAij, 



(45) 
(46) 



(d t ~C P )K = -D l D t a + a ( Aij A* 3 + -K z ) , (47) 



(fit - Cp) A^ = e'^ (—DiBjOL + aRij) 



TF 



a 



„( 



R-AaA 13 



KAa - 2A a A\ 



-K' 



(48) 



Notice that a=l will correspond to the case when the 
Ricci scalar is eliminated from the evolution equation for 
A^. 

As before, we will now concentrate on the case of 
geodesic slicing a = 1 with zero shift /3 l = 0, and consider 
a linear perturbation of flat space. 



g%j &ij ~t~ hij : 

The evolution equations then become 

d t ^ = -if/6, 
dthij = 2 Aij, 

d t K = 0, 
d t Aij = i?g } - SyRW (1 - a) /3, 

with Ri] the linear Ricci tensor. Notice now that to 
linear order K does not evolve at all: to linear order 



(49) 

(50) 
(51) 
(52) 
(53) 



the evolution of the gauge variables {</>, K} is therefore 
completely trivial. In particular, if K is chosen to be 
zero initially, it will remain exactly zero: no need for any 
exact cancellation. 

Now, quite generally the Ricci tensor can be separated 
into 



Rij = Rij 



Rl 



(54) 



The first term above Rij is the Ricci tensor associated 
with the conformal metric which to linear order is 



R 



(i) 



1/2 (Vg 



at 



(55) 



with the Ti defined just as before, but now using the 
conformal metric 



fi := d k hik-l/2dih. 



(56) 



The second term in (|34|) is the part of the Ricci tensor 
coming from the conformal factor <j) which to first order 
is 



R 



■tj '' - ~ 2 {(htijV ~ "'I S Uat' 

Notice that one can easily prove that 
det gij = 1 =>■ h = 0, 



(57) 



(58) 



so we could in principle eliminate the second term in 
Eq. d56|). As we will see below, this is a bad idea, so 
here we will just add instead a parameter £ that will be 
equal to if we eliminate h, and equal to 1 if we don't 
(but see the next section, where the T's are promoted 
to independent variables). We can then rewrite the first 
order Ricci tensor as 



R 



(i) 



-1/2 
- 2 [didj 



(59) 



Using this we can find the linearized version of the 
constraints 



where now 



djfi = 0; (hamiltonian) 

i 

dtfi = 0. (momentum) 
fi := djhjj - 8di(f). 



(60) 
(61) 

(62) 



As before, having found the linearized form of the evo- 
lution equations, we will proceed to make a Fourier analy- 
sis of the system. We then assume that we have a solution 
of the form 



6 



i(ujt — kx) 



i{ujt — kx) 

i 

K = K e l{ut - kx \ 

A.. — A., p i{u>t-kx) 

The evolution equations for <fi and hij imply 



K = — 6i uic 



(63) 
(64) 
(65) 
(66) 



(67) 



Substituting this into the evolution equations for K 

(68) 



and Aij we find 



u 2 h = k 2 Mh, 
where now h is a seven-dimensional vector 



0) hxxi hyyi h Z zi hxy-> faxz-i hy 



(69) 



and 





( 











\ 






TO22 


^23 


m 24 







m 3 i 


TO32 


^33 


m 34 





M = 


m 4 i 


TO42 


m 43 


m 44 









































V o 











1/ 


mix 


= 8- 


16(1 - 


-a)/3, 




TO22 ; 


= (£- 


1)(1 


-(1- 


-*)/3), 


TO 2 3 : 


= ^-24 


= e- 


-U + 


1)(1 


-a)/3, 


TO31 = 


= m 4 i 


= 4- 


- 16(1 


-a)/3, 


™32 : 


= ^42 




K-i 


(1- 


a)/3, 


™33 : 


= 777,44 


= 1 - 


-(€ + 


1)(1 


- CT )/3, 


m 34 = 


= ^43 




+ i 


(1- 


a)/3. 



(70) 



with 

(71) 
(72) 
(73) 
(74) 
(75) 
(76) 
(77) 

The hamiltonian and momentum constraints now re- 
duce to the three equations (again, not four!) 



h xx - 80 = 0, 


(78) 




(79) 


h xz = 0, 


(80) 



where as before, the first condition results from both the 
hamiltonian constraint and the x component of the mo- 
mentum constraints. 

The eigenvalues A and eigenvectors of the matrix (|70| ) 
are now somewhat more complicated. Let us consider 
first the eigenvalues on their own. They are: 

• A = 0, with multiplicity 3. 

• A = 1, with multiplicity 2. 



A = (a - 1 + 3£er ± rj) /6, with 



l + o- (34 - 420 + a 2 (1 + 3C)" 



1/2 



(81) 



There are a couple of things to notice from the last 
two eigenvalues. First, notice that if we take cr=0, one 
of these eigenvalues is always negative, which implies the 
existence of an exponentially growing mode, i.e. we have 
an unstable system of equations. So we must add some 
multiple of the hamiltonian constraint to the evolution 
equation of A^. How much we need to add will depend 
on the value of £. Moreover, with a little algebra one 
can also see that taking £=0 results as well in a nega- 
tive eigenvalue. This means that if we had decided to 
use the constraint h=0 (£=0) in the expression for the 
Ricci tensor, we would again have an unstable system of 
evolution equations. A safe value for £ turns out to be 
£=1. If we choose this, the characteristic structure of the 
matrix (|7(]) becomes 



(82) 
(83) 
(84) 
(85) 



(86) 
(87) 



• A = 0, with corresponding eigenvectors 

vi = (1,8,-4,-4,0,0,0), 
v 2 = (0,0,0,0,1,0,0), 
v 3 = (0,0, 0,0, 0,1,0), 
v 4 = (0,1,0,0,0,0,0). 

• A = 1, with eigenvectors 

v 5 = (0,0,1,-1,0,0,0), 
v 6 = (0,0,0,0,0,0,1). 

• A = (4a — l)/3, with eigenvector 

v 7 = (0, (4a + 2), (4a - 1), (4a - 1), 0, 0, 0) . 



Notice the last eigenvalue A = (4cr— l)/3 will only 
be positive for <r > 1/4, which tells us that we must 
add at least this much of the hamiltonian constraint to 
the evolution equation for Aij. A natural choice is to 
take er=l. This corresponds to completely eliminating 
the Ricci scalar from this equation, and results in the 
eigenvalue reducing to the physical speed of light. 

The type of solutions that the different eigenvectors 
represent are: 

1. Two physical solutions that travel with the speed 
of light (A = 1) represented by the transverse- 
traceless vectors V5 and vq. 

2. One mode that violates the hamiltonian constraint, 
the x component of the momentum constraint, and 
the constraint h = 0, that travels with the speed 
equal to the square root of (4er — l)/3, represented 
by the vector v-j. 
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3. Two modes that violate only the momentum con- 
straints, and "travel" with speed zero (A = 0) rep- 
resented by the vectors v-i and V3. 

4. One mode that violates the hamiltonian constraint, 
the x component of the momentum constraint, and 
the constraint h = that has speed zero (A = 0) 
represented by the vector V4. 

5. One pure gauge mode (satisfying all the con- 
straints) that travels with speed zero (A = 0) rep- 
resented by the vector v±. 

The structure of these eigenvalues and eigenvectors 
tells us in the first place that one has to be careful in 
the way in which different constraints are added to the 
evolution equations. The simple statement that one is in 
principle free to add multiples of constraints to evolution 
equations is true only if one does not worry about the 
stability of the hnal system. In this case we have seen 
how using blindly the constraint h—0 to simplify one of 
the equations results in the appearance of an unstable 
mode, and how neglecting to use the hamiltonian con- 
straint in another equation also gives rise to an unstable 
mode. A similar point has also been made in |4Cf| in the 
context of adding multiples of the hamiltonian constraint 
to the standard ADM evolution equations. 

From the characteristic structure described above, we 
can see that we now have four zero speed modes instead of 
three (assuming we do take £=1), so the situation would 
seem worse than before. Three of these modes are con- 
straint violating, and we will worry about them in the 
next section. What about the gauge mode? The gauge 
mode is of course still there, and it still has zero speed (as 
it should), but now it is in a much more convenient form. 
From looking at v\ we see that its evolution depends on 
the evolution equation for </>, which we have seen is triv- 
ial in the linear and non-linear case, and the evolution 
of the traceless part of hij , which is also trivial as long 
as the constraint trA=0 is satisfied (see Eq. ([[l])). The 
important point is the following: the fact that this mode 
evolves trivially is now the consequence of the simple al- 
gebraic constraint trA=0, and is independent of exact 
cancellations in derivatives of the metric that appear in 
the Ricci tensor. This provides an easy way to control the 
mode: Numerically setting trA to zero after each step of 
the evolution ensures that the gauge mode cannot grow. 

A comment is in order here. It has been recognized 
for some time [ pi)| , [il"lJl§| ] that gauge modes can propa- 
gate with arbitrary speeds. The analysis presented above 
shows that constraint violating modes can do the same. 
Often one does not think about these modes because they 
are unphysical, and one can avoid to excite them with 
an appropriate choice of initial data. However, from a 
numerical point of view, they will never really vanish 
and as we have just seen they can have important con- 
sequences on the stability of our evolutions. Even when 



these modes have a real speed of propagation (as opposed 
to an imaginary speed indicating an instability on the an- 
alytic level), if that speed is larger than the speed of light 
they can cause numerical instabilities if one forgets about 
their existence and chooses a time step based only on the 
extension of the physical light-cones. 

V. DEALING WITH THE CONSTRAINT 
VIOLATING MODES: USING THE MOMENTUM 
CONSTRAINTS 

In the previous section we have seen how separating 
out the gauge variables {0, K} provides a way to control 
the zero speed gauge mode. This still leaves us with the 
zero speed constraint violating modes to worry about. 
Here we will show how those modes can be dealt with by 
using the momentum constraints to modify the evolution 
equations of extra first order degrees of freedom. 

The idea of using the momentum constraints to modify 
the evolution equations is at the core of many recent hy- 
perbolic reformulations of the Einstein equations [p~p| p~3|| . 
In particular, the use of the momentum constraints to ob- 
tain evolution equations for extra first order variables can 
be traced back to the Bona-Masso formulation Jl0| , pl| . 
Here we will follow for simplicity the approach of Baum- 
garte and Shapiro p9|] (a very similar approach has been 
used before by Shibata and Nakamura p8|]). 

We will again concentrate on the case of geodesic slic- 
ing a = 1 with zero shift f3 l — 0, and consider a lin- 
ear perturbation of flat space. The linearized evolution 
equations were given by (p>0|)-(|53])). The Ricci tensor that 
appears in the evolution equation for Ay was separated 
as 

(^) 

with 

Rf) = -1/2 (VLA, - - djti) , (90) 

and 

4W=-2(W + %VLt^)- (91) 

Now, instead of writing the quantities Ti in terms of 
their definition ( |56| ) as we did before, we will promote 
them to independent quantities, and use ( p6| ) only to 
obtain their initial values. We will then need an evo- 
lution equation for the Ti. This we can obtain trivially 
from @ and (|l|): 

k 

Notice that we can use the fact that Ay is supposed to 
be traceless to eliminate the last term above. However, 



8 



we still don't know if this will turn out to be a good idea 
or not, so instead we again introduce a parameter £ and 
write 



(93) 



There is still one extra modification we want to make 
to this evolution equation: We will add to it a multiple 
of the momentum constraints (61) to obtain 



dtti = 2 (to - 1) d kA k i 



_ 4-777 

fflitrA - — 8iK, (94) 



with to arbitrary. Equation ( p4|) above is our final evolu- 
tion equation for the IV Keeping the I\ as independent 
variables, we also have to remember that we have intro- 
duced the extra constraints Ti = J2k ^fc^ifc- 

For the Fourier analysis, we again consider plane waves 
moving along the x direction. From the evolution equa- 
tions for 



and hij we find 



K = — 6i u/<f>, 



Aiq UjlljAA. 



(95) 



Substituting this in the evolution equations for Ti we 
obtain 



(m - 1 + <e/2) h xx 

- 8m(f> 



T x = —ik 

+ [hyy + h 

t y = -ik (to - 1) h xy , 
t z = -ik (to - 1) h xz . 



(96) 

(97) 
(98) 



And finally, substituting all these results back into the 
evolution equations for K and Aij we find 

w 2 h = k 2 Mh, 



where h is the same as before 
h 



hxx? hyy? h ZZl h X y, h xz , hy Z 

and where the matrix M is now 

/ 000 0\ 

TO21 TO22 m 2 3 TO24 

TO31 TO32 m-33 TO34 

M = TO4I TO42 ?7T-43 TO44 



(99) 
(100) 



TO 
TO 

1/ 



(101) 



with 



TO21 
m 2 2 
m 2 3 
m 3 i 
m 32 

TO34 



8 (1 - 2m) - 16 (1 - to) (1 - cr)/3, 
(2m + C-l)(l-(l-<r)/3), 
m a4 = e-(C + l)(l-ff)/3 > 
m 4i = 4 - 16 (1 - to) (1 - cr)/3, 
m 42 = -(2m + e-l)(l-fT)/3, 
m 44 = 1 - (£ + 1) (1 - tr)/3, 
TO43 = - (f + 1) (1 - ff)/3. 



(102) 
(103) 
(104) 
(105) 
(106) 
(107) 
(108) 



Notice that introducing the f j as independent vari- 
ables by itself does not change our analysis based on M, 
which is obtained by eliminating the fV But the evo- 
lution equations for the Ti motivate the introduction of 
the parameter to, whose effect we consider now. The 



eigenvalues of the matrix (101) turn out to be 



• A = 0, with multiplicity 1. 

• A = to, with multiplicity 2 

• A = 1, with multiplicity 2. 

• A = (1/6) \b± (b 2 -c) 1/2 



with 



b = -1 + a + 3£cr + 2to (2 + a) 
c = 36ct(-1 + C + 2to). 



(109) 
(110) 



The last two eigenvalues are quite complicated, so we 
will concentrate for the moment on the particular case 
(7=1. In that case the eigenvalues and eigenvectors of M 
reduce to 

• A = 0, with corresponding eigenvector 

vi = (1,8,-4,-4,0,0,0). (Ill) 

• A = to, with corresponding eigenvectors 

v 2 = (0,0,0,0,1,0,0), (112) 
v 3 = (0,0,0,0,0,1,0). (113) 

• A = 1, with eigenvectors 

v 4 = (0, 2£/(2 - 2to - 0, 1, 1, 0, 0, 0) , (114) 
« 5 = (0,0, 1,-1, 0,0,0), (115) 
v 6 = (0,0, 0,0, 0,0,1). (116) 

• A = 2to + £ — 1, with eigenvector 

v 7 = (0,1, 0,0, 0,0,0). (117) 

And the type of solutions represented are: 

1. Two physical solutions that travel with the speed 
of light (A = 1) represented by the transverse- 
traceless vectors v$ and ^6- 

2. One mode that violates the hamiltonian constraint, 
the x component of the momentum constraints, and 
the constraint h = that also travels with the speed 
of light (A = 1) represented by the vector w 4 . 

3. Two modes that violate only the momentum con- 
straints, and travel with speed to 1 / 2 represented by 
the vectors v 2 and v 3 . 
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4. One mode that violates the hamiltonian constraint, 
the x component of the momentum constraints and 
the constraint h = that has speed (2m + £ — l) 1 / 2 
represented by the vector vj. 

5. One pure gauge mode (satisfying all the con- 
straints) that travels with speed zero (A = 0) rep- 
resented by the vector V\. 

Notice first that all constraint violating modes have 
now acquired a non-zero speed. If we want to have all 
eigenvalues non- negative (and hence all speeds real), we 
must have 

m > 0, (118) 

and 

2m + £ - 1 > =>. m > i_ £. (119) 

In particular, if we take £ =0 (that is if we use the fact 
that trA=0 in the evolution equation for 1^) then we 
must have m > 1/2. So in order to have a stable sys- 
tem we must add a finite multiple of the momentum 
constraints to the evolution equation for I\. If we fail 
to use the momentum constraints, the system will have 
an exponentially growing mode. This is consistent with 
the results of the last section, where we didn't have the 
Fj (which in some sense is equivalent to not using the 
momentum constraints), and we found that taking £=0 
resulted in an unstable system. 
Notice also that if we take 

m = 1, £ = 0, (120) 

then we have one zero speed mode and six modes that 
travel with the speed of light. This is precisely the choice 
made by Baumgarte and Shapiro in p9| , so the result 
above explains why it was necessary in their case to add 
a multiple of the momentum constraints, and also why 
one should expect to have only the speed of light as a 
characteristic speed in their system. In the case m=l, 
£=0, the eigenvector v± might appear at first sight to be 
singular, but from the form that the matrix M takes in 
this particular case it is not difficult to show that in fact 
V4 is replaced by (0, 0, 1, 1, 0, 0, 0) with all other eigenvec- 
tors remaining unchanged. The only zero speed mode left 
is the pure gauge mode v\, but as we have seen before, its 
evolution does not rely any more on exact cancellations 
in the Ricci tensor. 

Finally, let us consider again the case when a ^ 1 , but 
now keeping m=l and £=0. In this case the eigenvalues 
and eigenvectors become 

• A = 0, with corresponding eigenvector 

wi = (1,8,-4,-4,0,0,0). (121) 



• A = 1, with corresponding eigenvectors 



v 2 = (0,0,0,0,1,0,0), (122) 

v 3 = (0,0,0,0,0,1,0), (123) 

v A = (0,-2, 1,1, 0,0,0), (124) 

v 5 = (0,0, 1,-1, 0,0,0), (125) 

v 6 = (0,0, 0,0, 0,0,1). (126) 

• A = a - , with eigenvector 

v? = (0,1, 1,1, 0,0,0). (127) 



We see now that depending on how large a multiple of 
the hamiltonian constraint we add to the evolution equa- 
tion of Aij, we can change the speed of propagation of 
the mode that represents the trace of hij (and hence the 
trace of Aij). If we do not use the hamiltonian constraint 
at all (cr=0), we will again have a zero speed unphysical 
mode. However, this is not as bad as it might seem be- 
cause in practice this mode is very easy to control since it 
will vanish if one imposes the algebraic constraint trj4=0. 

VI. NUMERICAL EXAMPLES: STABILITY OF 
MINKOWSKI SPACETIME 

To compare the stability properties of the different sys- 
tems in a simple situation we will consider the evolution 
of Minkowski spacetime, with a flat initial slice, but with 
a non-trivial spatial coordinate system. Since the extrin- 
sic curvature is zero, the spacetime should then remain 
static. Numerically, of course, the Ricci tensor will not be 
exactly zero, so we can expect some non-trivial evolution, 
but if the system is stable we will only have spurious nu- 
merical noise that should propagate away. If the system 
is unstable, however, we can expect that the numerical 
noise will slowly grow in amplitude. We will be evolving 
the full non-linear equations, so the initially slow growth 
of the numerical noise can be expected to trigger non- 
linear growth at late times. 

In order to obtain our initial metric, we start from the 
flat space metric in spherical coordinates 

dl 2 = dr 2 + r 2 dn 2 , (128) 

with dVt 2 the solid angle element. We then make the 
following coordinate transformation 

r = f(l-af(f)) , (129) 

with < a < 1 and f(r) a smooth monotonously de- 
creasing function that is 1 for small f and for large f. 
The particular form of the function / that we will use 
here is a Gaussian 

/(f) = e- f2 /° 2 . (130) 
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In terms of the new radial coordinate the metric be- 
comes 



with 



dl = .911 dr + f <?22 dQ 



911 = \l-a(f + rf)}\ 



922 = (1 - afY 



(131) 

(132) 
(133) 



Finally, for our 3D evolutions we transform this metric 
to Cartesian coordinates in the standard way, 



x = r sin f cos < 
y = f sin 9 sin q 
z = fcosO. 

So our initial metric is 



(134) 
(135) 
(136) 



9 xx 


= [x 2 9n 


f (y 2 + z 2 ) 


922 


If 2 , 


(137) 


9yy 


= [y 2 gn ■ 


f {x 2 + z 2 ) 


922 


If 2 , 


(138) 


9zz 


= [z 2 9n ■ 


f (x 2 + y 2 ) 


922 


If 2 , 


(139) 


9xy 


= xy (.gii 


-922) /f 2 , 






(140) 


9xz 


= xz (g n 


- .922) If 2 , 






(141) 


9yz 


= yz (511 


- 922) If 2 . 






(142) 



We must also say something about the gauge condi- 
tions used. For simplicity, we will use a zero shift vec- 
tor. For the lapse we could try geodesic slicing, but even 
small numerical perturbations will cause focusing (we are 
evolving the full non-linear Einstein equations) . It is bet- 
ter to use a slicing that can react to the evolution and 
can propagate away spurious numerical waves. Harmonic 
slicing is ideal for our purposes. It is defined via the fol- 
lowing evolution equation for the lapse 



dt,a = 



-a 2 K. 



(143) 



Since K is initially set to 0, the lapse should remain 1 if 
the evolution is exact. 

Finally, a comment about boundary conditions. We 
have used a very simple 'zero order extrapolation' bound- 
ary condition, that is, we update the boundary by just 
copying the value of a given field from its value one grid 
point in (along the normal direction to the boundary). 
This condition is not very physical, nor does it allow 
waves to leave the computational grid cleanly enough, 
but it is very robust, and can be used with all the dif- 
ferent formulations studied here in a stable way (at least 
for the time scales under study). Since our emphasis is 
on the stability of the interior evolution, we are content 
with having a stable boundary condition. We have used 
more sophisticated boundary conditions in various cases, 
but it is difficult to find one that will remain stable for 
all the evolution systems considered. 




FIG. 3. Surface plot of g X x along the x axis as a function of 
time for the simulation using the standard ADM formulation. 



We now present results of simulations performed with 
the different systems. The numerical method used in 
all these simulations was the so-called 'iterative Crank- 
Nicholson' (ICN) scheme with three iterations. We have 
found that three iterations are enough to obtain a stable, 
second order accurate numerical scheme [p2| . 

First we show the results of a simulation using the stan- 
dard ADM formulation for the case when a=0.1 and a=2. 
For this simulation we used a grid with 64 3 points and a 
resolution on Ax = 0.2. Figure || shows a surface plot of 
g xx along the x axis as a function of time. We see that g xx 
keeps its initial shape for some time, but at late times it 
starts to fall apart near the center. The simulation finally 
crashes at £=79. Figure |J shows the root mean square 
(r.m.s.) of the hamiltonian constraint over the numerical 
grid as a function of time. We see that for a long time 
there is an essentially linear growth of the r.m.s. of the 
hamiltonian constraint superimposed with small oscilla- 
tions, just what we expect from the linear analysis of the 
previous sections. At late times, however, the non-linear 
effects take over and we have a catastrophic blow-up, as 
we argued above. 

Next, we show results of the conformally rescaled sys- 
tem of section [V, using £=1, and two different values 
of a: cr=0 (no use of the hamiltonian constraint) and 
(7=1 (use of the hamiltonian constraint to completely 



eliminate the Ricci scalar from the evolution equation 
for Aij). From our analysis we expect the system with 
cr=0 to have an exponentially growing mode and thus to 
be very unstable. The a— I should only have the zero 
speed modes and should be much more stable (but still 
not completely stable). Figure || shows the r.m.s. of the 
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FIG. 4. Root mean square of the hamiltonian constraint as 
a function of time for the simulation using the standard ADM 
formulation. 




FIG. 5. Root mean square of the hamiltonian constraint as 
a function of time for the simulation using the standard CT 
formulation with £=0 and two different values of a. 

hamiltonian constraint for these two runs. We see that 
our predictions are indeed correct, the tr=0 run becomes 
rapidly unstable and crashes at i=4, while the a=l is far 
more stable and only crashes at i=33. 

We now show the results of the choice <r=0, m=l, £=0 
in section |v], as used by Baumgarte-Shapiro pgj] . We 
have set trAij to zero at each step as discussed above. 
Figure |^ shows again a surface plot of g xx along the x 
axis as a function of time (but notice the change of scale) . 
The evolution now goes past i=500 with no trace of an 
instability. Figure [j] shows the r.m.s. of the hamilto- 
nian constraint for this run. The hamiltonian constraint 
rapidly becomes much larger than in the ADM case at 
early times (by almost a factor of 10). However, it then 
stops growing and simply oscillates around a constant 
value, showing again no sign of the linear growth or the 
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FIG. 6. Surface plot of g xx along the x axis as a function 
of time for the simulation using the Baumgarte-Shapiro for- 
mulation. 

blow-up that we saw for ADM. 

Finally, we show results of a series of simulations 
done by keeping <r=0 and £=0, but changing the value 
of m (the amount of momentum constraint added to 
the evolution equation of the F's). Figure ^ shows 
the r.m.s. of the hamiltonian constraint for runs with 
to = {0,0.25,0.5,0.75} (compare with the m=l case 
shown above). As expected from our analysis, we see 
that the cases with to < 1/2 rapidly become unstable. 
The simulation with to=0 crashes at f=4 while the one 
with to=0.25 crashes at £=12. On the other hand, the 
cases with to > 0.5 remain stable past i=400. 

VII. CONCLUSIONS 

We have studied the stability properties of the stan- 
dard ADM formulation of general relativity based on a 
linear perturbation analysis. We focus attention on the 
zero speed modes. We conjecture that the zero speed 
modes can cause instabilities in evolutions of the ADM 
system in its standard form. These instabilities do not 
have a numerical origin, but rather they correspond to 
genuine blowing-up solutions of the differential equations. 

We show that the zero speed modes come in two forms: 
a pure gauge mode that satisfies all the constraints, and 
is therefore a legitimate physical solution, and a series of 
non-physical constraint violating modes. We investigate 
the change in behavior of these modes going from the 
standard ADM formulation to the conformal-traceless 
(CT) systems of Shibata and Nakamura |2l| and Baum- 
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Norm of hamiltonian constraint 



FIG. 7. Root mean square of the hamiltonian constraint 
as a function of time for the simulation using the Baum- 
garte-Shapiro formulation. 



garte and Shapiro p9| , and their derivatives. Two fea- 
tures we believe responsible for the better stability prop- 
erty of the conformal systems are identified: 1. The zero 
speed gauge mode is governed by an equation that is free 
from the complication of the Ricci tensor, thus decou- 
pling it from the rest of the system. 2. The constraint 
violating zero speed modes, on the other hand, acquire 
a finite speed of propagation due to the introduction of 
extra first order degrees of freedom, and the use of the 
momentum constraints to modify the evolution equations 
for these degrees of freedom. We present numerical ex- 
amples to support our analysis. 

We consider the study presented in this paper as a first 
step towards the understanding of the stability issue in 
the numerical evolution of the Einstein equations. 
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APPENDIX A: FINITE DIFFERENCE 
APPROXIMATION TO THE LINEARIZED ADM 
EQUATIONS 



E 0.0020 




FIG. 8. Root mean square of the hamiltonian constraint 
as a function of time for the simulation using the Baum- 
garte-Shapiro formulation with different multiples of the mo- 
mentum constraint added to the evolution equation for the 
F's (different values of the parameter m). 



We will consider a simple finite difference approxima- 
tion to the linearized ADM evolution equations written 
in second order form. For this we start from equations (j^) 
and (Q), and substitute one into the other to find 



dlKj - Vfl at /iij = didjh -2^2 d(id m h 



(Al) 



We now construct a simple second order finite differ- 
ence approximation to this equation using standard cen- 
tered differences, 



1 



2 fii 



(Aty 
i 



s 2 f 

2 u t Jm ! 



u i J m i 



(A2) 
(A3) 



(A4) 
(A5) 



(Ax) 2 tJm 
with f m = f(t = nAt, Xi = mi Ax) and 

r2 rn fn+1 o rn , pn—1 

r2 rn rn o rn , rn 

°i Jm ~ Jnii + 1 L Jm,i T J mi -1 

Let us now consider a plane wave solution of the form 

(hij) m = h ij e i(nuAt+m - kAx ^ . (A6) 

But notice now that we allow the waves to move along 
any direction. This is because even if different directions 
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are equivalent from the analytic point of view, they are 
not equivalent numerically because the numerical grid 
introduces preferred directions. 

If we substitute this into the finite difference approxi- 
mation to (uu we find the following equation, 



cos(uAt)] h = Mh 

P 



(A7) 



where p := At/ Ax is the Courant parameter, h is defined 
as before, 



hxxi h"yy> h ZZl h X yj h X z> hyz I ; 



(A8) 



and M is the matrix 



M 



( 


u\ + u 2 z 


ul 




ISxy 


^>Sxz 









u 2 


u% + u 2 z 




ISxy 





1Sy Z 






K 


u\ 


Ul+U 2 y 





ISxz 


ISyz 












— S X y 


u\ 


Syz 


$xz 














Syz 


u l 






\ 


~ s yz 








$xz 


$xy 


U 2 X 


) 



where we have defined 

u\ := 2 [1 - cos(fc l Aa;)] , 
Sij := —sia(kiAx)sbx(kjAx) 

Let us now define 
2 



A 



[1 - cos(wAi)] 



Equation (A7) now becomes 
Mh = Ah 



(A9) 



(A10) 
(All) 



(A12) 



(A13) 



which is just an eigenvalue equation. Here we face one 
problem: the characteristic polynomial is of 6th order, 
and is difficult to solve exactly in the general case. We 
will then consider a couple of particular cases. 

First, assume that the wave moves only on the x direc- 



tion, so k y 



0. In this case everything simplifies 



considerably, and we find that the eigenvalues of M are 
just 

• A = 0, with multiplicity 3. 

• A = u 2 , with multiplicity 3. 

This has precisely the same structure we found before 
for the exact system of differential equations. The only 
difference being that the wave speed is now not quite 1. 
The wave speed in fact depends on the wave number k x , 
and can be obtained from the dispersion relation 



cos(kjAi)] 



(A14) 



Notice that for small k x Ax (large wavelengths com- 
pared to the grid spacing) this relation reduces to 



U3 



(A15) 



which is what one expects. For smaller wavelengths we 
obtain wave speeds that are smaller than 1, showing the 
dispersive nature of the finite difference approximation. 

The results above are not particularly surprising. One 
obtains essentially the same thing for the simple wave 
equation. The interesting case is when we consider waves 
moving in a direction different from the coordinate lines. 
We will then consider the particular case of waves moving 
in the diagonal direction, for which k x = k y — k z = 
k. The characteristic polynomial now does not simplify 
nearly as much, but one can still find the eigenvalues 
analytically. They are 



• A = 

• A = 



.A = | 



— s 2 , with multiplicity 2. 

+ 2s 2 , with multiplicity 2. 

5u 2 - 2s 2 ± (9u 4 - 12s 4 + 12w 2 s 2 ) 



1/2 



where 



2 [1 - cos(fcAa;)] 



s 2 = sin 2 (fcAa;). (A16) 



The values of the different roots are shown in Fig- 
ure |jj The solid lines indicate the four distinct eigen- 
values, while the dashed line indicates the eigenvalue one 
would obtain (also along the diagonal line) for the finite 
difference approximation to the simple 3D wave equa- 
tion A = 3u 2 . The plot is only in the region fcAx 6 [0, n] 
since larger wave numbers can not be represented on the 
numerical grid (k = tt/Ax is the so-called 'Nyquist' fre- 
quency of our grid). 

There are several things to notice from this result. 
First, we now have four distinct eigenvalues instead of 
two: the numerical grid has broken the degeneracy of 
the exact problem. Second, the three eigenvalues that 
where zero in the exact case are now only zero for k = 0, 
and are clearly non-zero for any finite k. This shows that 
the zero speed modes have picked up a non-zero speed in 
the numerical approximation. This artificial speed is very 
small for large wavelengths (small fc), but becomes con- 
siderable for smaller wavelengths. Finally, we see that for 
small values of k we recover the exact result: one eigen- 
value vanishes as fc 6 /12, two as fc 4 /4, and the other three 
go to zero as 3fc 2 , which is the correct result for waves 
traveling with a speed of 1 along the diagonal. 
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